Visualize Cell Health Predictions on Repurposing Hub Data

Gregory Way, 2019

suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(ggplot2))

source(file.path("scripts", "visualize_utils.R"))

Set Constants

consensus <- "modz"

output_dir <- file.path("figures", "umap", consensus)
dir.create(output_dir)

Warning message in dir.create(output_dir):
“'figures/umap/modz' already exists”

Load and Process Data

# Load Data
data_dir <- file.path("repurposing_cellhealth_shiny", "data")
real_file <- file.path(
    data_dir,
    paste0("moa_cell_health_", consensus, ".tsv.gz")
)

cp_embedding_df <- readr::read_tsv(real_file, col_types = readr::cols())

cp_embedding_df <- cp_embedding_df %>%
    dplyr::mutate(Metadata_Treatment = "Compound")

cp_embedding_df$Metadata_Treatment[cp_embedding_df$Metadata_broad_core_id == "DMSO"] = "DMSO"

print(dim(cp_embedding_df))
head(cp_embedding_df, 3)

[1] 10368    88

A tibble: 3 × 88

Metadata_Plate_Map_Name	Metadata_broad_core_id	Metadata_broad_sample	Metadata_dose_recode	Metadata_mmoles_per_liter	Metadata_pert_well	umap_x	umap_y	broad_id	pert_iname	⋯	cell_health_modz_target_cc_g2_high_h2ax	cell_health_modz_target_cc_late_mitosis_n_spots_h2ax_mean	cell_health_modz_target_cc_cc_high_h2ax	cell_health_modz_target_vb_percent_dead_only	cell_health_modz_target_cc_s_high_h2ax	cell_health_modz_target_cc_cc_n_spots_h2ax_per_nucleus_area_mean	cell_health_modz_target_cc_s_n_spots_h2ax_per_nucleus_area_mean	cell_health_modz_target_cc_g1_plus_g2_count	cell_health_modz_target_vb_live_cell_width_length	Metadata_Treatment
<chr>	<chr>	<chr>	<int>	<dbl>	<chr>	<dbl>	<dbl>	<chr>	<chr>	⋯	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<chr>
C-7161-01-LM6-001	BRD-A25234499	BRD-A25234499-001-18-3	6	10.0000	B13	-1.580029	-0.9219316	BRD-A25234499	aminoglutethimide	⋯	0.2751957	0.2319843	0.03473731	0.2267767	0.237375	0.10548474	0.1536775	-0.5309856	0.14445485	Compound
C-7161-01-LM6-001	BRD-A25234499	BRD-A25234499-001-18-3	5	3.3333	B14	-1.904766	-2.4360995	BRD-A25234499	aminoglutethimide	⋯	0.5169502	0.2421589	0.01257728	0.6328903	0.399858	0.02291539	0.2027969	-0.1597692	-0.29454513	Compound
C-7161-01-LM6-001	BRD-A25234499	BRD-A25234499-001-18-3	4	1.1111	B15	-1.899164	-2.5914881	BRD-A25234499	aminoglutethimide	⋯	0.4012496	0.2682545	0.08088620	0.2074322	0.337209	0.18061697	0.2859738	-0.2541855	-0.04831094	Compound

Visualize Metadata and Cell Health Variables

ggplot(cp_embedding_df,
       aes(x = umap_x, y = umap_y)) +
    geom_point(aes(color = Metadata_dose_recode,
                   size = paste(Metadata_Treatment)),
               pch = 16,
               alpha = 0.6) +
    theme_bw() +
    scale_color_viridis_c(name = "Dose\nRecoded") +
    scale_size_discrete("Treatment",
                        range = c(0.5, 2)) +
    xlab("UMAP 1") +
    ylab("UMAP 2")

output_file <- file.path(
    output_dir,
    paste0("umap_repurposing_cell_painting_dose_consensus_", consensus, ".png")
)
ggsave(output_file, height = 5, width = 6, dpi = 500)

Warning message:
“Using size for a discrete variable is not advised.”

ggplot(cp_embedding_df,
       aes(x = umap_x, y = umap_y)) +
    geom_point(aes(color = Metadata_dose_recode,
                   size = paste(Metadata_Treatment)),
               pch = 16,
               alpha = 0.6) +
    theme_bw() +
    scale_color_viridis_c(name = "Dose\nRecoded") +
    scale_size_discrete("Treatment",
                        range = c(0.5, 2)) +
    xlab("UMAP 1") +
    ylab("UMAP 2")

output_file <- file.path(
    output_dir,
    paste0("umap_repurposing_cell_painting_dose_consensus_", consensus, ".png")
)
ggsave(output_file, height = 5, width = 6, dpi = 500)

Warning message:
“Using size for a discrete variable is not advised.”

ggplot(cp_embedding_df %>% dplyr::filter(Metadata_Treatment == "DMSO"),
       aes(x = umap_x, y = umap_y)) +
    geom_point(aes(color = Metadata_pert_well),
               pch = 16,
               size = 2,
               alpha = 0.6) +
    geom_point(data = cp_embedding_df,
               color = "grey",
               alpha = 0.1,
               size = 0.5) +
    theme_bw() +
    theme(legend.position = "none") +
    xlab("UMAP 1") +
    ylab("UMAP 2")

output_file <- file.path(
    output_dir,
    paste0("umap_repurposing_cell_painting_dose_consensus_dmso_", consensus, ".png")
)

ggsave(output_file, height = 5, width = 6, dpi = 500)

Certain Models had Better Performance in A549

We applied all models to the Drug Repurposing Set data.

Here, output visualizations of several of the top models.

# Load feature mapping
mapping_dir <- file.path("..", "1.generate-profiles", "data", "labels")
mapping_file <- file.path(mapping_dir, "feature_mapping_annotated.csv")
map_df <- readr::read_csv(
    mapping_file,
    col_types = readr::cols(.default = readr::col_character())
)

print(dim(map_df))
tail(map_df, 3)

[1] 75 16

A tibble: 3 × 16

id	readable_name	original_name	feature_type	measurement	gate_required	assay	hoechst	edu	ph3	gh2ax	caspase	draq7	cell_rox	dpc	description
<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>
cc_g2_g1_count	Cell Count - G2/G1	G2/G1	cell_cycle	cell_cycle_count	1	hoechst_edu_ph3	1	1	1	0	0	0	0	0	G2/G1
cc_g1_plus_g2_count	Cell Count - G1+G2	G1+G2	cell_cycle	cell_cycle_count	1	hoechst_edu_ph3	1	1	1	0	0	0	0	0	G1+G2
cc_g2_plus_all_m_count	Cell Count - G2 + M	G2 + All M-phase	cell_cycle	cell_cycle_count	1	hoechst_edu_ph3	1	1	1	0	0	0	0	0	G2+M

Visualizing Specific Cell Health Models

Only a few are shown here, but all are saved in a separate folder.

Predicting Number of Live Cells

visualize_umap(
    cp_embedding_df,
    target_variable = "cell_health_modz_target_vb_num_live_cells",
    legend_title = "Num Live Cells",
    output_dir = "none",
    save = FALSE,
    print_figure = FALSE
)

Predicting Live Width:Length

visualize_umap(
    cp_embedding_df,
    target_variable = "cell_health_modz_target_vb_live_cell_width_length",
    legend_title = "Live Cell\n(Width:Length)",
    output_dir = "none",
    save = FALSE,
    print_figure = FALSE
)

Predicting Live Cell Roundness

visualize_umap(
    cp_embedding_df,
    target_variable = "cell_health_modz_target_vb_live_cell_roundness",
    legend_title = "Live Cell Roundness",
    output_dir = "none",
    save = FALSE,
    print_figure = FALSE
)

Predicting Number of Objects

visualize_umap(
    cp_embedding_df,
    target_variable = "cell_health_modz_target_cc_all_n_objects",
    legend_title = "Number of Objects",
    output_dir = "none",
    save = FALSE,
    print_figure = FALSE
)

Predicting Live Cell Area

visualize_umap(
    cp_embedding_df,
    target_variable = "cell_health_modz_target_vb_live_cell_area",
    legend_title = "Live Cell Area",
    output_dir = "none",
    save = FALSE,
    print_figure = FALSE
)

Predicting Number of Objects in Cell Cycle Stage

visualize_umap(
    cp_embedding_df,
    target_variable = "cell_health_modz_target_cc_cc_n_objects",
    legend_title = "Num Cell\nCycle Objects",
    output_dir = "none",
    save = FALSE,
    print_figure = FALSE
)

Predicting Number of Objects in G1

visualize_umap(
    cp_embedding_df,
    target_variable = "cell_health_modz_target_cc_g1_n_objects",
    legend_title = "G1 Objects",
    output_dir = "none",
    save = FALSE
)

Predicting EdU positive intensity

visualize_umap(
    cp_embedding_df,
    target_variable = "cell_health_modz_target_cc_s_intensity_nucleus_area_mean",
    legend_title = "Sum S phase",
    output_dir = "none",
    save = FALSE
)

Predicting ROS

visualize_umap(
    cp_embedding_df,
    target_variable = "cell_health_modz_target_vb_ros_back_mean",
    legend_title = "ROS Background",
    output_dir = "none",
    save = FALSE
)

Visualize All Cell Health Variables

cell_health_variables <- colnames(
    cp_embedding_df %>%
        dplyr::select(starts_with("cell_health_modz_target_"))
    )

length(cell_health_variables)

70

pdf_file <- file.path(
    output_dir,
    paste0("repurposing_hub_umaps_consensus_", consensus, ".pdf")
)
pdf(pdf_file, width = 5, height = 5, onefile = TRUE)

for (cell_health_variable in cell_health_variables) {
    umap_gg <- visualize_umap(
        df = cp_embedding_df,
        target_variable = cell_health_variable,
        legend_title = "Prediction:",
        title = cell_health_variable,
        dpi = 200,
        save_figure = FALSE
    )
}

dev.off()

pdf: 2

Visualize Model Scores on DMSO samples

assay_theme_file <- file.path("..", "3.train", "scripts", "assay_themes.R")
source(assay_theme_file)

col_types <- readr::cols(
    .default = readr::col_character(),
    shuffle_false = readr::col_double(),
    shuffle_true = readr::col_double()
)

rank_file <- file.path(
    "repurposing_cellhealth_shiny",
    "data",
    paste0("A549_ranked_models_regression_", consensus, ".tsv")
)
model_rank_df <- readr::read_tsv(rank_file, col_types = col_types)

# Recode the target variable
model_rank_df$target <- paste0("cell_health_", consensus, "_target_", model_rank_df$target)

head(model_rank_df, 3)

A tibble: 3 × 10

metric	target	original_name	readable_name	feature_type	measurement	assay	description	shuffle_false	shuffle_true
<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<dbl>	<dbl>
r_two	cell_health_modz_target_vb_live_cell_width_length	Live Width:Length	Live Width / Length	viability	shape	dpc	Width/Length	0.8916483	6.635413e-03
r_two	cell_health_modz_target_vb_num_live_cells	# Live Cells	# Live Cells	viability	death	draq	Number of live cells	0.8811216	-5.934792e-06
r_two	cell_health_modz_target_vb_live_cell_roundness	Live Cell Roundness	Live Cell Roundness	viability	shape	dpc	Cell Roundness	0.8552069	2.012126e-01

dmso_embeddings_df <- cp_embedding_df %>%
    dplyr::filter(Metadata_Treatment == "DMSO")

non_dmso_embeddings_df <- cp_embedding_df %>%
    dplyr::filter(Metadata_Treatment != "DMSO")

std_dev_dmso_features <- apply(
    dmso_embeddings_df %>% 
        dplyr::select(matches("cell_health_modz_target")),
    2,
    sd
)
std_dev_compound_features <- apply(
    non_dmso_embeddings_df %>%
        dplyr::select(matches("cell_health_modz_target")),
    2, 
    sd
)

std_dev_all_df <- dplyr::bind_cols(
    as.data.frame(std_dev_dmso_features),
    as.data.frame(std_dev_compound_features)
) %>%
    dplyr::mutate(
        features = colnames(dmso_embeddings_df %>%
                                dplyr::select(matches("cell_health_modz_target")))
    ) %>%
    dplyr::left_join(model_rank_df, by = c("features" = "target")) 

good_performing <- std_dev_all_df %>%
    dplyr::filter(shuffle_false > 0)

bad_performing <- std_dev_all_df %>%
    dplyr::filter(shuffle_false <= 0)

std_dev_good_df <- good_performing %>%
    dplyr::mutate(performance_scaled = (
        good_performing$shuffle_false - min(good_performing$shuffle_false)
    ) / (
        max(good_performing$shuffle_false) - min(good_performing$shuffle_false)
    )
                  )

print(dim(std_dev_good_df))
head(std_dev_good_df, 2)

[1] 37 13

A data.frame: 2 × 13

	std_dev_dmso_features	std_dev_compound_features	features	metric	original_name	readable_name	feature_type	measurement	assay	description	shuffle_false	shuffle_true	performance_scaled
	<dbl>	<dbl>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<dbl>	<dbl>	<dbl>
1	0.05967878	0.2168907	cell_health_modz_target_cc_late_mitosis_n_spots_h2ax_per_nucleus_area_mean	r_two	pH3 neg, Hoechst cond (late mitosis) - Number of Spots per Area of Nucleus - ...	Late M - # of gH2AX Spots per Area of Nucleus	cell_cycle	dna_damage	hoechst_edu_ph3_gh2ax	In late M cells: average number of gH2Ax spots per nucleus area	0.01027311	-0.08045787	0.001896058
2	0.35065337	1.5849058	cell_health_modz_target_cc_g1_n_objects	r_two	G1 - Number of Objects	G1 - # cells	cell_cycle	g1_phase	hoechst_edu_ph3	Number of G1 cells	0.77424668	0.06178493	0.867049811

ggplot(std_dev_good_df,
       aes(x = std_dev_dmso_features, y = std_dev_compound_features)) +
    geom_point(aes(color = assay, size = performance_scaled),
               alpha = 0.8) +
    geom_point(data = bad_performing,
               aes(color = assay),
               size = 1,
               alpha = 0.5) +
    theme_bw() +
    geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "red") +
    scale_color_manual(name = "Assay",
                       values = dye_colors,
                       labels = dye_labels) +
    scale_size_continuous(name = "Performance Scaled") + 
    xlab("Standard Deviation across DMSO Well Consensus") +
    ylab("Standard Deviation across all Consensus Compound Treatments") +
    theme(axis.text = element_text(size = 8),
          axis.title = element_text(size = 10),
          strip.text = element_text(size = 6),
          legend.text = element_text(size = 6),
          legend.title = element_text(size = 8),
         legend.key.size = unit(0.4, "cm"))

output_file <- file.path(output_dir, "dmso_vs_compound_standard_deviation.png")
ggsave(output_file, height = 5, width = 6, dpi = 500)